## Honest DID Analysis (Results are estimated by lfe package version 3.0-0)
## These results are reported in Figure 5
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate long and short sample
data5=subset(data4, T>40)
data_base1=data4
data_base1$D=0
data_base1$D[which(data_base1$Group=="Treated")]=1
data_base2=data5
data_base2$D=0
data_base2$D[which(data_base2$Group=="Treated")]=1

###########################
# Standard DID
twfe_results_AOD=fixest::feols(AOD ~ i(weeks, D, ref = 9)  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. 
                               | Year + weekday + weeks*Plant_Code + T,   cluster = "Plant_Code",
                               data = data_base1)
betahat <- summary(twfe_results_AOD)$coefficients[1:13] #save the coefficients
sigma <- summary(twfe_results_AOD)$cov.scaled[1:13,1:13] #save the covariance matrix


###########################
# Sensitivity analysis using relative magnitudes restrictions
delta_rm_results <-
  HonestDiD::createSensitivityResults_relativeMagnitudes(
    betahat = betahat, #coefficients
    sigma = sigma, #covariance matrix
    numPrePeriods = 8, #num. of pre-treatment coefs
    numPostPeriods = 5, #num. of post-treatment coefs
    Mbarvec = seq(0.1,0.6,by=0.1) #values of Mbar
  )
delta_rm_results
originalResults <- HonestDiD::constructOriginalCS(betahat = betahat,
                                                  sigma = sigma,
                                                  numPrePeriods = 8,
                                                  numPostPeriods = 5)

HonestDiD::createSensitivityPlot_relativeMagnitudes(delta_rm_results, originalResults)

